# This R code was build under version 4.1.1 
# R packages needed
# plm (version 2.4.3) for lineal panel data models 
library(plm)
# stargazer (version 5.2.2) for tables
library(stargazer)
# sandwich (version 3.0.1) for robust standard errors
library(sandwich)

#Files "panel.csv" and "cross.csv" must be in the same path as this file.

# reading panel data
panel <- read.csv("panel.csv")

#table 2a: panel summary statistics

stargazer(panel[,-c(1,2,3,5,12)], type="html", 
          omit.summary.stat=c("p25", "p75"), 
          title="Table 2a: Dependent and control variables descriptive statistics", 
          covariate.labels = c("Area", "Population", "Conflict Intensity", "Coca cultivated area","Per-capita tax revenue",
                               "IDPs expulsion rate", "Collective titles granted", 
                               "Land suitability", "Land suitability dummy"),
          out="table2a.doc" )



# reading cross-section data
cross <- read.csv("cross.csv")

#table 2b: cross-section statistics
stargazer(cross[,c(12,13,4)], type="html", 
          omit.summary.stat=c("p25", "p75"), 
          title="Table 2b: Dependent and control variables descriptive statistics", 
          covariate.labels = c("Productive land under collective tenure (2014)", "Agricultural land area (2014)"),
          out="table2b.doc" )

#table 2 = table2a+table2b


#function to compute robust standard errors
rse <- function(nomb){
  covv <<- vcovHC(nomb, type="HC1")
  assign(paste0('rse_',gsub(pattern = "r_", replacement = "", substitute(nomb))),sqrt(diag(covv)), envir = .GlobalEnv)
}

#basic panel regressions (table 3)

r_all1 <- plm(Exp_rate~ lag(collect_land), 
                data= panel,
                index=c("MUNICIP_code","year"), 
                effect ="twoway",
                model = "within")

r_all2 <- plm(Exp_rate~ lag(collect_land)+pop+conflict_int+Coca_hect+tax_rev_per_mill, 
                data= panel,
                index=c("MUNICIP_code","year"), 
                effect ="twoway",
                model = "within")

r_all3 <- plm(Exp_rate~ lag(collect_land)*Land_suitab+pop+conflict_int+Coca_hect+tax_rev_per_mill,

                data= panel,
                index=c("MUNICIP_code","year"), 
                effect ="twoway",
                model = "within")

r_all4 <- plm(Exp_rate~ lag(collect_land)*Land_suitab_dummy+pop+conflict_int+Coca_hect+tax_rev_per_mill,
                
              data= panel,
              index=c("MUNICIP_code","year"), 
              effect ="twoway",
              model = "within")

# output for table 3
stargazer(r_all1, r_all2, r_all3, r_all4, type = "html", out="table3.doc",
          se = list(rse(r_all1), rse(r_all2), rse(r_all3), rse(r_all4)),
          covariate.labels = c("Collective titles granted(t-1)", "Population", "Conflict Intensity", "Coca cultivated area",
                               "Per-capita tax revenue", "Collective titles granted(t-1)*Land suitability", 
                               "Collective titles granted(t-1)*Land suitability dummy"),
          keep.stat = c("n","rsq")
          )

# panel regressions robustness checks (table 4)

r_allr1 <- plm(Exp_rate~ lag(collect_land)*Land_suitab+lag(pop)+lag(conflict_int)+lag(Coca_hect)+lag(tax_rev_per_mill),
         
               data= panel,
               index=c("MUNICIP_code","year"), 
               effect ="twoway",
               model = "within")

r_allr2 <- plm(Exp_rate~ collect_land*Land_suitab+pop+conflict_int+Coca_hect+tax_rev_per_mill, 
                
                data = panel,
               index=c("MUNICIP_code","year"), 
               effect ="twoway",
               model = "within")

r_allr3 <- plm(Exp_rate~ collect_land+pop+conflict_int+Coca_hect+tax_rev_per_mill, 
                 
                 data = subset(panel, pacifico>0),
               index=c("MUNICIP_code","year"), 
               effect ="twoway",
               model = "within")

# output for table 4
stargazer(r_allr1, r_allr2, r_allr3, type = "html", out="table4.doc",
          keep = c("collect_land","lag(collect_land):Land_suitab","collect_land:Land_suitab"),
          se = list(rse(r_allr1), rse(r_allr2), rse(r_allr3)),
          covariate.labels = c("Collective titles granted(t-1)", "Collective titles granted(t-1)*Land suitability", 
                               "Collective titles granted",
                               "Collective titles granted*Land suitability"),
          keep.stat = c("n","rsq")
          )
#cross-section regressions (table 5)


r_cross_1<- lm(Exp_rate_2015~prod_area_col_title2014, data=cross)

r_cross_2<-lm(Exp_rate_2015~prod_area_col_title2014+pop2014+conflict_int+Coca_hect+tax_rev_per_mill_2014+agric_area_2014+Area+Altitude, data=cross)
r_cross_3<- lm(Exp_rate_2015~prod_area_col_title2014+pop2014+conflict_int+Coca_hect+tax_rev_per_mill_2014+agric_area_2014+Area+Altitude 
               +Atlantica +Oriental+Pacifica+Central, data=cross)
r_cross_4<- lm(Exp_rate_2014~prod_area_col_title2014+pop2014+conflict_int+Coca_hect+tax_rev_per_mill_2014+agric_area_2014+Area+Altitude 
               +Atlantica +Oriental+Pacifica+Central, data=cross)

# output for table 5
stargazer(r_cross_1, r_cross_2, r_cross_3, r_cross_4, type = "html", out="table5.doc",
          se = list(rse(r_cross_1), rse(r_cross_2),rse(r_cross_3), rse(r_cross_4)),
          covariate.labels = c("Productive land under collective tenure 2014","Population", "Conflict Intensity",
                               "Coca cultivated area","Per-capita tax revenue 2014","Agricultural land area 2014","Area",
                               "Altitude"
                               ),
          keep.stat = c("n","rsq")
          )



